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Abstract 


Models currently used for analyses of thermal and water behavior of a PEM fuel cell are based 3D computational fluid dynamics (CFD). 
However, the analyses are limited to a single cell with static behavior. Thus, these models cannot be used for analyses of dynamic behavior of a 
stack that continuously varies according to operating conditions. The model proposed describes dynamic behavior of a stack with two adjoining 
cells and endplate assembly, and work as a current controlled voltage source that can be used for optimization of BOPs and the associated controls. 
Simulations have been conducted to analyze start-up behaviors and the performance of the stack. Our analyses deliver following results: (1) dynamic 
temperature distribution in both the through-plane direction and the along channel direction of the fuel cell stack, (2) effects influencing the source 
terms of current density, and (3) dynamic oxygen concentration distribution. The temperature profile and its variation propensity are comparable 
to the previous results [Y. Shan, S.Y. Choe, J. Power Sources, 145 (1) (2005) 30-39; Y. Shan, S.Y. Choe, J. Power Sources, in press]. 


© 2006 Elsevier B.V. All rights reserved. 
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1. Introduction 


Current computational models available in both the academic 
world and the market are either too simple or complex to particu- 
larly describe dynamic behaviors of a PEM fuel cell stack. Some 
authors [1—3] simply employ empirical equations, whose coeffi- 
cients are obtained by fitting a polarization curve. This approach 
is useful for a design of the power system, but ignores effects 
of temperature, water and reactants on the cell performances. 
Thus, it can hardly describe the complex behaviors of a stack. In 
addition, 2D or 3D models using CFD techniques proposed by 
other authors [4—7] can capture the complexity of a single cell, 
but are limited to steady state analyses and unable to represent 
an unsteady behavior of a stack. Um et al. [8] published a 2D 
CFD model that describes the transient behavior of the bulk flow, 
species and electro-chemical reactions in a single cell, which has 
been extended to an isothermal 3D model [10,11]. Likewise, 
Dutta et al. [9] developed a 3D model with an isothermal flow in 
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a cell that embeds a serpentine-type gas channel. However, the 
simulation is conveyed by the use of commercial software pack- 
age, Fluent, and they studied both the gas distribution and water 
generation in a single cell. Um and Wang [10] developed a 3D 
CFD model for a single cell. They intensively studied species and 
water removal in a straight and an inter-digitated flow channel 
and found enhancement of the performance at the inter-digitated 
shape. Furthermore, Wang and Wang [11] simulated a single cell 
with 36 gas serpentine channels taking a low humidity condi- 
tion by using the software package of Star-CD and presented the 
mechanism of the species transport and the associated current 
density distribution. Unlike the studies above, Ju et al. [12] con- 
sidered a non-isothermal condition in the 3D plane and simulated 
a cell with a straight channel by using Star-CD. 

All of works above, however, have been focused on descrip- 
tion of a single cell and are still unable to describe a stack 
and time-varying behaviors. On the other hand, the dynamic 
behavior of a stack can be improved by adding a simplified ther- 
modynamic model, which is proposed by Sundaresan [13,14]. 
The model regards a cell as a composition of layers and is used 
to analyze the start-up behavior from a sub-freezing tempera- 
ture. However, the model does not fully consider several factors: 
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Nomenclature 


membrane extension coefficient (m? m^?) 
mole concentration (mol m^?) 

thermal capacity (J mol! K^!) 
diffusion coefficient (m? s7!) 

Faraday number (C mol!) 

current density (A m2) 

current density (A m7?) 

mass flow rate (kg m?) 

thermal conductivity (W m`! K7!) 

gas permeability (m7) 

length (m) 

mole mass (kg mol!) 

electro-osmotic coefficient 

universal gas constant (J mol`! K7!) 
source term 

temperature (K) 

velocity (m s7!) 

voltage (V) 

dimensionless mol concentration (m?) 


qe Gas Si cR Ee 


Greek letters 


a transfer coefficient 

y water transfer coefficient (m s7!) 
E porosity 

n over-potential (V) 

À water content 

u viscosity (kg m^! s7!) 
p density (kg m^?) 

o conductivity (s m^ 1) 
[o potential (V) 
Subscripts 

a anode 

c cathode 

e electrolyte 

OC open circuit 

S solid matrix 

tota anode total 

totc cathode total 
Superscripts 

eff effective 

ref reference value 


(1) flow rate of species at the inlet of the channel must be the 
same as that at the outlet of the channel. Thus, no fluid dynam- 
ics are considered; (2) heat source terms in both the catalysts 
are empirically calculated with values suggested by the Wohr 
and Peinecke's model [15]. Accordingly, the anode source term 
is presumed as a relatively large value that in fact should be 
referred to be around zero [16]. As a result, the model does not 
show asymmetric phenomena of performance through the stack. 
Wetton et al. [17] proposed an explicit stack thermal model 
with the coolant channel coupled with a 1D cell model [18]. 


It shows a great temperature gradient of the stack, but with no 
dynamics at all. We proposed an enhanced quasi 1D stack model 
[19,20] based on the previous single cell model that considers 
the thermal and fluid dynamics. As a result, the model proposed 
is capable of capturing the dynamic temperature distribution 
including the asymmetrical effects in the stack, but missing the 
water distribution that are improved by adding an empirical rela- 
tionship between the flooding effect and the current density and 
temperature. 

As the matter of fact, none of current models can fully 
describe the stack behavior. On the system aspects, the model 
should be a current controlled voltage source, so that the load can 
be easily integrated into a stack model. Moreover, the domain 
should be so set up that can integrate two endplates and bus 
plate at the anodic and cathode side with an interface plate that 
embeds a coolant channel as well as two bipolar plates with the 
repeating basic cell unit. 

In fact, the mass and charge transport as well as the heat flux 
in the basic cell unit is described by the use of the Navier-Stokes 
and the potential and energy conservation equation. Likewise, 
the heat flux in the coolant channel is described by using both the 
Navier-Stokes and the energy conservation equations, while the 
rest of plate regions are described by the heat conservation equa- 
tions. The partial differential equations (PDEs) are solved by the 
SIMPLE [21] algorithm. In the following sections, details on the 
2D models are summarized, which includes assumptions, sim- 
ulation set-up and equations used for a description of the model 
proposed. Included are the procedure to solve the equations and 
a generation of grids. At the end, simulations are conveyed with 
boundary conditions used for the individual domain and the 
results are discussed. 


2. Model description 
2.1. Modeling domain and assumptions 


Major assumptions are made for a 2D stack model as follows: 


ja 


. Reactants as ideal gases; 

2. Incompressible and laminar flow; 

3. Isotropic and homogeneous electrodes, catalyst layers and 

membrane; 

4. Identical inlet conditions of each cell for both the cathode 
and anode as well as coolant channel; 

. Constant thermal conductivity of the materials in a fuel cell; 

. Neglect diffusions caused by multi-component; 

. No contact resistance; 

. No liquid water generated; 

. The source/sink term can be neglected in a PEFC where 
electrochemical reactions occur [26]. 


NO 00 -1 Q0 tA 


In addition, it is assumed that a single cell has a structure 
of sandwiched layers shown in Fig. 1. The anode sides of the 
cells are located on the left hand side, while the cathode sides 
on the right hand side. The single cell domain for the model 
is constructed with seven different layers that are symmetrically 
placed at the membrane layer. A gas flow channel, a gas diffusion 
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Fig. 1. Single cell schematic configuration [8]. 


layer and a catalyst layer for the anode side are located at the 
left side of the membrane layer as well as those for cathode side 
with a reversed order. Thus, a stack can be easily constructed by 
repeating this basic unit domain and adding bipolar plates with 
coolant channel, interface and bus plates, and end plates shown 
in Fig. 2. 

Finally, a stack model is completed by coupling of the 
domains for the basic units with two endplate assemble. Finally, 
simulation can be performed. 


2.2. Model description 


2.2.1. Charge transport 

Protons and electrons are the positively and negatively 
charged ions. The proton transfer in the proton conducting 
regions and the electron transfer in the electronic conducting 
regions determines the potential distribution in a cell. The poly- 
mer as electrolyte in the membrane and the catalysts belongs 
to the proton conducting region, while the catalysts, GDLs and 
BPs including gas flow and coolant channels are regarded as 
electrode, which is repeating in the stack configuration. 


Endplate 


Bus and interface (IF) 


Coolant channel 


Gas flow channel (oxidant) 
Diffusion layer 

Catalyst layer (Cathode) 
Membrane 

Catalyst layer (Anode) 


Diffusion layer 


Gas flow channel (Fuel) 


Coolant channel 


Bus and interface 


Endplate 


Fig. 2. Stack schematic configuration. 


The potentials in the electrolyte are governed by the potential 
conservation equation according to the Ohm’s law: 


V(ae V Pe) + So — 0 (1) 
So = jinthetwo catalyst layers (1a) 
So = Oin the membrane layer (1b) 


where the proton conductivity oe is a function of temperature 
and the water content in the polymer material [24]; 


1 1 
= 100 1268 | —— — 
m s ( e 3) 


x (0.005139AH,0/so;, — 0.00326) (2) 


According to the Butler- Volmer equation [1,8], the current 
densities in the anode and cathode catalysts can be expressed by 
the exchange current density, reactant concentration, tempera- 
ture and over-potentials according to the Butler- Volmer Eqs. (3) 
and (4): 


: j Xy 1/2 fa, +a 
ja = aif ( 2 ) ———— Fn (3) 
Xp, ref RT 
x 1/2 Oc F 
; ref. Oo c 
=— € exp| — 4 
Jc ai (22) »( i) (4) 


where the surface over-potential is defined as a difference 
between the electrodes and electrolyte referring to an equilib- 
rium state. 


n(x, y)=V— Vequlibrium = ®; — d. — Voc (5) 


where ®, and ®e are the potentials of the electrons conducting 
solid materials and electrolyte, respectively, at electrodes and 
electrolyte interface. The open circuit potential at the anode is 
assumed to be zero, while the open circuit potential at the cathode 
becomes a function of a temperature as [1,9]: 


Vic = 0.00257 + 0.2329 (6) 


Then, the local current density for the protons can be simplified 
with 


I = —o.V®_ (7) 


On the other hand, the electronic conducting region includes 
the entire stack except the membranes and the two endplates. In 
fact, the electronic conductivity in the entire electron conducting 
layers is at least two orders higher than the proton conductiv- 
ity. Thus, the potential drop caused by the electrons transfer is 
negligible. Therefore, the values of the potentials in the electron 
conducting regions can be regarded as a single value. When the 
load current is applied to the model as an input, the voltage drops 
in the electron conducting regions vary because of the flow field 
being changed. This effect is reflected by using a current con- 
servation equation at an equilibrium state of the potential field 
and described by the following equations: 


f(D) = f i, dV — TinputLch = Oin anode catalyst layers (8) 
V 
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f(E) = f i, dV — TinputLch = Oin cathode catalyst layers 
V 
(9) 


Accordingly, both of the electrolyte potential distribution and 
electrode potential (®, and d) can be corrected. The resulting 
equations are implicit and can be solved numerically. 


F(®,. + ABe) = i, Ja (de+4A8.) dV — linpuLcn = 0 (10) 
V 


F(®, + A@s) = | Jc(@s+.A6,) dV — linputLcn = 0 (11) 
v 


Then, the new electrolyte potential is obtained by adding the 
correction factor calculated in the anode catalyst layer to the 
previous value of the potential for the domain of the catalyst and 
the membrane, while the previous of the potentials are obtained 
by solving the potential conservation equation (Eq. (1)). Like- 
wise, the new electrode potentials are calculated by adding the 
correction factor calculated in the cathode catalyst layer to the 
previous value of the potential. 


PEW = G+ Ade, PY = pM t Ad, ad 


2.2.2. Mass transport 

2.2.2.1. Anode/cathode side. The anode side includes a gas 
channel, a gas diffusion layer, and catalysts. Like the anodic 
side, the cathode side includes a cathode gas channel, a gas dif- 
fusion layer and a catalyst. The role of the layers is to provide 
a pathway for the transfer of the fuel and humidity from the 
inlet to the catalysts. The flow field can be a shape of either 
serpentine or inter-digitated channels. Compared to the serpen- 
tine shape, the inter-digitated type can dramatically increase the 
mass transfer in the gas diffusion layer by convection. In this 
study, a straight channel is selected to represent a serpentine 
channel. Then, the mass transport in the channel is governed by 
the following equations: 


e Mass conservation 


ð 3 
(05) | voci) = 0 (13) 
ot 
e Momentum conservation 
(peù ER = 
px + V(peliii) = —eV p + View" vii) + Su (14) 
Su = ral in the gas diffusion and catalyst layers [5] 
(14a) 
Su = 0 in the gas channel (14b) 


Similarly, the hydrogen concentration can be obtained by using 

the species conservation equation: 

O(eX Hp ) 
ot 


+ V(eŭ Xn) = V(Dip, V Xu) + Sm (15) 


SH = — B in gas diffusion layers and catalyst layers [1] 
2 FCtota 
(15a) 
SH, =0 in gas channel (15b) 
where the effective diffusivity Dil is, 
Di - e^ Dy, Xo, = m (16) 
Ctota 
and the species conservation equations for cathode are: 
o(eX 3 
ae + V(eiXo,) = VDEVX0,) + So, (17) 
Je 
So. = in catalyst layers (17a) 
? 4 FCtotc i A 
So, =0 in gas channel and gas diffusion layers (17b) 
9(eXnjo) 2 
— y + VGRXno)- V(ODigoVXio)- Smo — (18) 
Jc À 
Smo = — in catalyst layers (18a) 
3 2 F'Ctotc 
Smo —O in gas channel and gas diffusion layers (18b) 


where the effective diffusivities and mole fractions ofthe oxygen 
and water are as follows: 


ff 1.5 eff 1.5 CO» 
Do, = £ ` Do, Dio —£"Dmo, Xo, = —, 
Ctotc 
CH30 
Xpmo- — (19) 
Ctotc 


2.2.2.2. Water transport in membrane. Water distribution in the 
membrane is determined by the electro-osmotic force, the dif- 
fusion and convection. At a single phase, the influence of the 
convection is negligible and the water concentration can be 
described by the following equations: 


H20 
0c. ? 
ot 


I 
= V(D{O vce?) — y (a (20) 


The first term describes the diffusion dependent water flux and 
the second one does the water transport dependent upon the 
electro-osmotic force. And the electro-osmotic coefficient is a 
function of the local water content (nq): 


ÀH50/S05 
22 


The water content Ag, 0/so; is a function of the water uptake 


obtained from the water concentration: 
H20 
Ce 


oy [M° — bibo 


na=2.5 (21) 


ÀH20/S0; = (22) 

The water flux at the boundary interfaces between the cat- 
alysts and membrane is described by using the Robin type 
boundary condition, which presents the relationship between 
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the water concentration in the gas phase and in the membrane 
[25]. 


c e e,eq 


JEO = yO _ 0) _ T4 Q3) 


However, compared to the reactants transport phenomena in 
the gas diffusion layer, the water transport in the membrane 
is more complicated. Some experiments showed the different 
time scales of the hydration and dehydration processes in the 
Nafion [25]. It takes about tens of seconds for liquid water to 
get adsorbed, which is two orders of magnitude larger than for 
the dehydration process. This hydration-dehydration dichotomy 
might result from two reasons: (1) the water diffusion coeffi- 
cient in the membrane is a function of local water concentration. 
And it varies with the changing of water concentration; (2) the 
Robin type water transfer coefficient y is also affected by the 
water concentration at the boundaries between the membrane 
and catalysts. 

Springer et al. [23] proposed a coefficient of water diffusion 
equation that is derived from experimental results: 


1 1 
pH20 — p' 2416 | — — — 24 
e exp ( (xs =)) vm 


where 


2.64227E — 13X5,0/so, 
7.75E — 11Ag,0/so, — 9.5E — 11 
2.5625E — 11A5,0/so, — 2.1625E — 10 


D' = 


According to the data proposed Berg [7], the Robin type water 
transfer coefficient is chosen with yı 25 x 107^ when a mem- 
brane is fully immersed into liquid water, while yg =4.5 x 1076 
for the membrane dried in a gaseous condition. 


2.2.2.3. Coolant flows. The coolant flow channels embedded in 


Cell 1 


Endplate 
— — Busplate 
——— |Fplate 


—— — |Fplate 
— —- Busplate 
— —- Endplate 


0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 
X Axis 


Fig. 3. Stack geometry information. 


and the heat source has occurred at the wall, the coolant releases 
most of the heat to the environment via the radiator. If the coolant 


1.23 > ÀH20/S03 
6> ÀH20/S0; = 1.23 


(24a) 


14 > Amoj/so, > 6 


flow channel is assumed as a straight channel, the flow of the 
coolant can be governed by the following equations: 


e Mass conservation 


bipolar plates or endplates provide part of the pathway for the dp =) 

o = + V(ou)-0 (25) 
coolant circuits. After the heat exchange between the coolants ot 
Table 1 
Parameters for 2D models 
Quantity Value 
Gas channel length, L (cm) 7.112 
Oxygen diffusivity in gas (cm? s7!) 5.2197 x 107? 
Hydrogen diffusivity in gas (cm? s^!) 2.63 x 107? 
Dissolved oxygen diffusivity in active layer and membrane (cm? s7!) 2.0 x 1074 
Dissolved hydrogen diffusivity in active layer and membrane (cm? s7!) 2.59 x 1076 
Faraday constant, F (C mol!) 96,487 
Permeability of backing layer, K (cm?) 1.76 x 10-6 
Universal gas constant, R (J mol-! K7!) 8.314 
Cathodic transfer coefficient 2 
Anodic transfer coefficient 2 
Inlet nitrogen-oxygen mole fraction 0.79/0.21 
Air-side inlet pressure/fuel-side inlet pressure (atm) 5/3 
O» stoichiometric flow ratio 3.0 
Hp stoichiometric flow ratio 2.8 
Reference exchange current density x area of anode (A cm?) 5.0 x 107 
Reference exchange current density x area of cathode (A cm?) 1.0 x 1074 
Total mole concentration at the anode side (mol cm?) 66.81 x 1076 
Total mole concentration at the cathode side (mol cm?) 17.81 x 1076 
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Table 2 
Geometry parameters for the fuel cell 


Thickness, m 


Heat conductivity (W m^! K7!) 


Heat capacity (JIkg—! K7!) Density (kg m^?) 


GDL 0.0004 4 840 2000 
Catalyst layer 0.000065 0.2 TIO 387 
Membrane layer 0.000183 0.21 1100 1967 
Gas channel 0.001 52 935 1400 
Plate 0.001 52 935 1400 
Coolant channel 0.001 30 935 1400 
GDL porosity 0.4 
Catalyst layer porosity 0.2 
GDL tortuosity 3.725 
Bipolar plate contact area percentage 0.5 
Membrane molecular mass (kg mol! ) 1.1 
Fuel cell area (m?) 0.0367 
Fuel cell active area (m2) 0.03 
e Momentum conservation Table 3 
m Initial values 
O(pii) "m E 
or + V(puu) = —V p + V(uVu) (26) Quantity Value 
Temperature (K) 353 
2.2.3. Heat flux in a stack Oxygen nondimensinal concentration 0.21 
The heat in the stack is produced by five different sources Cathode inlet velocity (m s^!) 0.334 
Anode inlet velocity (m s7!) 0.157 


that include the entropy and losses caused by over-potentials at 
two catalysts, proton conductivity, electron conductivity, and the 
phase change of water. Under the condition of a single phase, 
there is no heat generation associated with the phase change. 
In addition, all the heat generated is then transferred from the 
source to the fluid by conduction and convection and completely 
removed out of the stack by the coolant and gases at the channel 
outlet. Then, the thermal behavior of a stack can be governed by 
the energy conservation equation [5]: 


9 EN > > 

3; 0m Cp T) + V(EpriCp T) = VRE VT) + Sr (27) 
. dVoc P. 

Sr=j(n+T + — inthe catalyst layers (27a) 

dT Oe 

P 

$T = in the membrane layer (27b) 
Oe 

Sr —0 inthe gas channel and gas diffusion layers (27c) 


where the overall density and thermal conductivity are defined 
as 


K — ke + gem Q8) 


stack 


Pm = Pf + Pstack: 


and the fluid mixture properties are 


pf = S aX: prCp,r = a Q9) 
i i 


2.2.4. Boundary conditions 

2.2.4.1. For the fluid flow. The fluid velocities at the inlet of the 
anode, the cathode, and the coolant channel are predetermined. 
The standard exit boundary and no-slip boundary conditions 
are applied to the channel exits and walls, respectively. In the 
species field, the inlet species concentrations are given, while 


the species gradients at the channel exits and walls are set to 
zero. All boundary conditions for the fluid flow are summarized 
as below. 

At the anode inlet (30) and outlet (31) 


u=0, v= Vp, in XH, = XM, ins TH, = Th; in (30) 
ð ax ita, 
“= 0, =i. Ses Ih _ 9 
dy dy dy 
Ty, 
m l0 (31) 
dy 


Pressure drop in a single cell at 1.0 sec 


p 


200000 
200000 
200000 
199999 
199999 
199999 
199999 
199999 
199999 
199998 
199997 
199996 
199996 
199996 
199996 
199995 
199995 


y (m) 


0.0615 0.062 0.0625 0.063 0.0635 
x (m) 


Fig. 4. Pressure drop in the cell 1 (Unit: Pa). 


202 
At the cathode inlet (32) and outlet (33) 
u=0, 


v= vo, Xo, = Xo,is Tair = 


air,in 
ð ox ox 
SU: LT: Cc IU cdd % o, 
dy dy dy 
OT air —0 
ay 


At the coolant inlet (34) and outlet (35) 


u =Q, U = UCoolant,in: TCoolant = TCoolant,in 
Qv aTe 
u=0, LR 0, coolant ih 
dy dy 
At the wall 
ox dXo, 
u = 0, v= 0, n > %2 ecl 
dy dy 
ox 
H20 =0 
dy 


(32) 


(33) 


(34) 


(35) 


(36) 


2.2.4.2. For the electrolyte potential field. The boundary con- 
ditions for the gradient of the electrolyte potential field are set 


Oxygen concentration at 0.3 sec 


y (m) 


0.056 


0.058 


0.06 
x (m) 


0.062 0.064 


Oxygen concentration at 1.0 sec 


0.056 


0.058 0.06 
x (m) 


0.062 0.064 
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to zero at the left anode catalyst as well as the right cathode 
catalysts 


(37) 


2.2.4.3. For the electrolyte potential field. The boundary con- 
ditions for the gradient of the electrical potential field is set as 
zeros at the right anode catalyst boundary as well as the left 
cathode catalyst boundary 


0d, 
ox 


=0 (38) 


3. Numerical solution 


First of all, the conservation equations are discretized by 
using the control-volume-based finite difference method. The 
flow solution procedure is based upon the SIMPLE algorithm 
[21] with a collocated grid cell centered approach. 

The simulation set-up for the stack includes two endplate 
assemblies, two coolant and flow channels, and two cells with 
one bipolar plate (Fig. 3). 

The input parameters used for the current simulation are 
summarized (see Tables 1-3). 


Oxygen concentration at 0.5 sec 


0.07 


gas2 
0.205 
0.195 
0.185 
0.175 
0.165 
0.155 
0.145 
0.135 
0.125 
0.115 
0.105 
0095 
0.085 
0.075 
0.065 
0.055 
0.045 
0.035 
0.025 
0.015 
0005 


0.056 0.058 0.06 


x (m) 


0.062 0.064 


Oxygen concentration at 10.0 sec 


gas2 
0.205 
0195 
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0.165 
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0.145 
0135 
—0.04 9.125 
0115 
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0095 
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x (m) 
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Fig. 5. Dynamics of oxygen concentration at two cells. 
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4. Analyses 


Fig. 4 shows the simulated result of the pressure field in the 
cell number 1. The pressure drop on the anode flow channel 
is very small and negligible, simply because the viscosity of 
the hydrogen is smaller than that of the air on the cathode. In 
addition, the velocity of the anode gas is slower than that of the 
cathode one. Itis noted that the Reynolds number for the cathode 
is 23.7 for the simulation and the resulting pressure drop is 10 Pa 
approximately, which is comparable to the results in Ref. [22]. 

Fig. 5 shows transient behaviors of the oxygen concentrations 
in the cathode side flow channels. When the simulation starts, 
the oxygen concentrations are changing, but still identical in the 


Cathode Catalysts current density cell 1 at 0.3 sec 


0.07 
cdc 
-5E+06 
0.06 -1E+07 
-1.5E*07 
-2E+07 
-25E*07 
0.05 -3E*07 
-3.5E+07 
-4E+07 
-4.5E*07 
€ 0.04 -5E«07 
> 
0.03 
0.02 
0.01 
0.05628 0.05629 0.0563 
x (m) 
Cathode Catalysts current density cell 1 at 1.0 sec 
0.07 
cde 
-5E+06 
0.06 -1E+07 
-1.5E*07 
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both channels because the temperature effects on the reaction 
is not so high and the reactants consumed are not so different 
in the two cells. However, the concentration on the left side is 
lower than the right side in the channel because of the oxygen 
being consumed by the chemical reactions. The steady state has 
been reached after the 1 s. 

Fig. 6 shows a transient behavior of the source term of the 
current density generated in the cathode catalysts of the cell 1 
at the 0.3, 0.5, 1 and 10s. The magnitude of the current den- 
sity largely varies in the 1s. In fact, the source term of the 
current density are generally influenced by three phenomena 
with the rise of the temperature, the distribution of the cath- 
ode overpotentials and the concentration of the reactants. The 
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Fig. 6. Current density on the cathode side (Unit: A m7’). 


204 


y (m) 


y (m) 


y (m) 


Y. Shan et al. / Journal of Power Sources 165 (2007) 196—209 


Potential of electrolyte in cell 1 at 0.1 sec 
0.07 


0.06 


0.05 


0.04 


0.03 


0.02 


0.01 


0.0561 


0.0562 
x (m) 


pop 
-0.00493239 
-0.00904981 
-0.0145066 
-0.0254744 
-0.0360557 
-0.0445299 
-0.05 
-0.0588105 
-0.067378 
-0.0707184 
-0.0745611 
-0.0823546 
-0.084113 
-0.0854694 
-0.0864321 
-0.0921249 
-0.0944224 
-0.0975198 
-0.100531 
-0.101968 
-0.103077 
-0.103166 
-0.103452 


Potential of electrolyte in cell 1 at 0.5 sec 


0.0561 


0.0561 


0.0562 
x (m) 


0.0562 
x (m) 


0.056: 


0.056 


pop 


-0.00493239 
-0.00904981 
-0.0145066 
-0.0254744 
-0.0360557 
-0.0445299 
-0.05 
-0.0588105 
-0.067378 
-0.0707184 
-0.074561 1 
-0.0823546 
-0.084113 
-0.0854694 
-0.0864321 
-0.0921249 
-0.0944224 
-0.0975198 
-0.100531 
-0.101968 
-0.103077 
-0.103186 
-0.103452 


Potential of electrolyte in cell 1 at 3.0 sec 


pop 
-0.00493239 
-0.00904981 
-0.0145086 
-0.0254744 
-0.0360557 
-0.0445299 
-0.05 
-0.0588105 
-0.067378 
-0.0707184 
-0.0745611 
-0.0823546 
-0.084113 
-0.0854694 
-0.0864321 
-0.0921249 
-0.0944224 
-0.0975198 
-0.100531 
-0.101968 
-0.103077 
-0.103166 
-0.103452 


y (m) 


y (m) 


y (m) 


Potential of electrolyte in cell 1 at 0.3 sec 


0.07 


0.06 


0.05 


0.04 


0.03 


0.02 


0.01 


0.0561 


0.0562 


x (m) 


0.0563 


pop 
-0.00493239 
-0.00904981 
-0.0145066 
-0.0254744 
-0.0360557 
-0.0445299 
-0.05 
-0.0588105 
-0.067378 
-0.0707184 
-0.074561 1 
-0.0823546 
-0.084113 
-0.0854684 
-0.0864321 
-0.0821249 
-0.0844224 
-0.0975198 
-0.100531 
-0.101968 
-0.103077 
-0.103166 
-0.103452 


Potential of electrolyte in cell 1 at 1.0 sec 


0.07 


0.05 


0.04 


0.03 


0.02 


0.01 


0.0561 


0.0562 
x (m) 


0.0563 


pop 
-0.00493239 
-0.00904981 
-0.0145066 
-0.0254744 
-0.0380557 
-0.0445299 
-0.05 
-0.0588105 
-0.067378 
-0.0707184 
-0.0745611 
-0.0823546 
-0.084113 
-0.0854694 
-0.0864321 
-0.0921249 
-0.0944224 
-0.0975198 
-0.100531 
-0.101968 
-0.103077 
-0.103166 
-0.103452 


Potential of electrolyte in cell 1 at 10.0 sec 


0.07 


0.06 


0.05 


0.04 


0.03 


0.02 


0.01 


0.0561 


0.0562 


x (m) 


Fig. 7. Membrane potential dynamics at the cell 1 (Unit: V). 
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Fig. 8. Stack temperature 2D (Unit: K). 
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temperature in the first 10s is not so high and the influence is 
negligible. The overpotential on the left side of the cathode cat- 
alyst in the cell 1 is higher than that on the right side because 
the protons flow in the direction where the electrolyte potential 
decreases. 

The concentration of the reactant at the inlet side is higher 
than that on the outlet side. Therefore, the current density gets 
higher at the inlet side. However, the magnitude of the influ- 
ences through and along the planes depends upon, which one 
is dominant in an operating condition. For an example, when 
the concentration of the reactants decreases along channel, the 
overpotential gets increased because of the change of the mem- 
brane conductivity that is explained in details by Ju et al. [12]. 
As a matter of fact, water generated is being accumulated at 
the outlet side and hydration rate in the membrane becomes 
higher, while the inlet side relatively gets dehydrated because 
of less water accumulated. Likewise, the concentration of the 
reactants gradually decreases through the planes from the GDL 
to the membrane, while the overpotential increases. As a conse- 
quence, the source term of the current density depends upon an 
instant which one is dominantly influenced. 

It is observed that the current density along the channel direc- 
tion is dynamically varying within the | s. The current density at 
the outlet side decreases, while the one at the inlet side increases. 
This effect is caused by the change of the concentration of the 
reactants along the channel. 

The results presented in this paper are not able to represent 
some effects associated with water concentration, overpoten- 
tials, concentration of the reactants caused by two-phase flow. 
In addition, due to the limited calculation time of the codes 
developed, simulation has been run for only 10s and the steady 
state is not reached yet. 

Fig. 7 shows the simulated results of a phase potential field 
in the electrolyte of the cell 1 representing a catalyst and mem- 
brane. The potential field has reached to a steady state at the 
first second, where the potential drop in the membrane along the 
channel direction is not high. In fact, the gradient of the mem- 
brane conductivity along the gas channel direction is negligible 
because of the constraints on the single phase and requires a long 
computational time to see the effects. Conversely, the potential 
in the cathode catalysts shows a high gradient, because the chem- 
ical reaction at the inlet side is high than those on the outlet side. 
Subsequently, protons at the inlet side are more consumed and 
tend to flow toward the inlet direction. 

When the reactant concentration is large, the source term of 
current density is bigger. As a result, the potential of the elec- 
trolyte at the inlet gets smaller, the over-potential gets smaller. 
It decreases the source term of the current density, finally an 
equilibrium state is reached. 

Fig. 8 illustrates a transient behavior of the temperature distri- 
bution for a stack including two cells at eight different instants: 
0.1, 0.3, 0.5, 1, 3, 5, 8, and 10s. When the stack gets operated, 
the heat is generated and temperature rises. 

At the first instant, the most heat is generated in two cat- 
alyst layers when the reaction begins and at the same time 
the entropy change occurs. Thus, the temperature peak appears 
in the cathode catalyst layers of two cells instead of the 
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Fig. 9. Fluid temperature. 


membrane layer. In addition, the source term of the current 
density at the inlet region is higher than that of the outlet 
region, so the temperature along the channel direction gets 
decreased. 

At the following instants, the heat flux begins to diffuse and 
temperature in all the planes of the stack rises. 

Fig. 9 shows the temperature distribution of the gases and 
fluids in the cell 1 that includes two gas flow channels. The 
geometry of the cell is referred to Tables 2 and 3. The tem- 
perature profiles in both of channels are not symmetrical and 
the temperature on the cathode side is higher than the anode 
side because of the larger amount of heat generated in the cath- 
ode than in the anode. Due to the large heat transfer area of 
the GDL by the porosity, the temperature gradient between the 
catalyst and the GDL is relatively low, while the temperature 
drop at the interface between the GDL and the flow channel is 
high. 

Fig. 10 shows the dynamics of the temperature distribution. 
Firstly, the temperature dynamically changes at different instants 
when the times go by from 0.1, 0.3, 0.5, 1, 3, 5, 8, and 10s. It is 
observed that the temperature in the cells at the very beginning 
rapidly rises from an initial value of 353 K because of the heat 
generated in the catalysts by chemical reaction and in the mem- 
branes by ohmic losses. The shapes of the rising temperature 
in both of cells are identical. It is shown that the peak value in 
the catalyst layer of the cell 1 gradually becomes higher than 
the one in the cell 2. In fact, the endplate assembly used for 
the stack should be comparably thick because of the mechanical 
requirement for robustness and represents two large heat sinks. 
However, the distance between the layer with the heat source 
and the end assemblies in the stack is different in a typical stack 
construction and consequently heat transfer properties are not 
identical to each other. Finally, the temperature profile through 
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Fig. 10. Temperature distribution through the plane. 


the plane gets asymmetrical and the one on the left side cell 


becomes higher. 


Fig. 11 shows the temperature distribution in the cells 1 and 
2, which allows for a direct comparison of how the temperature 
profiles are dynamically changing at a different instant. It is 


207 


shown that the shapes of the temperature profile at the 0.1 s are 
the same. The endplate effect is to recognize at the 0.2 s, where 
the cell 1 shows a low temperature at the anode side. At the 0.5 
and 1 s, the gap between two cells gets larger. The temperature 
of the cell 1 on the cathode side is lower than that of the cell 2, 
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Fig. 11. Temperature distribution through the plane. 


while the temperature of the cell 1 on the anode side is higher ate. The temperature of the cell 1 becomes higher than the 
than the one of the cell 2. cell 2. 

At the | and 3s, the gaps on the both sides get larger and At the 8 and 10s, the behavior of the temperature profiles 
the peak values of temperature for both cells begin to devi- remains as the same as before. As a consequence, the perfor- 
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mance of the cell 1 becomes better than the cell 2. However, 
there are some constraints on the analyses conveyed because of 
ignored effects of two-phase flow. 


5. Conclusion 


In this work, a new dynamic 2D model for the PEM fuel 
cell stack considering the fluid and thermal characteristics is 
proposed. Emphases are placed on numerically solving the equa- 
tions and effects of the temperature distribution on the stack 
performance, which varies dynamically during operations. The 
model developed assumes a structure sandwiched by layers that 
include membrane, catalysts, gas diffusion layers, bipolar plates 
and endplate assembly. The domains specified for the model- 
ing are determined by the way of physically working principles 
rather than the layers used commonly. The separate setup of 
modeling domains provided an easiest way to understand the 
physics involved and consequently reduce the development time 
for the codes. 

The model can calculate dynamic distribution of pressure 
and reactants, current density, temperature and potentials in a 
stack. Due to the consideration of two adjoining cells, the mod- 
els can be used for particularly designing a stack and BOPs 
along with controls, where the interacting influences from the 
neighboring cells plays a significant role. For example, calcula- 
tions of temperature distribution across the stack can contribute 
to develop a thermal management strategy for the coolant 
circuit. 

The analyses and modeling proposed lay a groundwork, 
which ultimately should increase efficiency and performance of 
the system. Future work will include an expansion of the model 
for two-phase flow in a 3D geometry and optimize the model 
with 1D and 2D. 
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